using DataAnalyticsTools.Models;
using System;
using System.Collections.Generic;
using System.Linq;

namespace DataAnalyticsTools.Core.ChangePointDetect
{
    /// <summary>
    /// 基于回归的变化点检测器（使用Chow Test）
    /// 通过Chow Test的F统计量检测回归结构突变，适合检测趋势斜率的变化
    /// </summary>
    public static class RegressionDetector
    {
        /// <summary>
        /// 最大置信度上限
        /// </summary>
        private const float MAX_CONFIDENCE = 0.95f;

        /// <summary>
        /// 使用Chow Test检测变化点
        /// </summary>
        /// <param name="values">待检测的时间序列数据</param>
        /// <param name="sensitivityFactor">灵敏度因子（默认2.25，建议范围0.9-6.0）
        /// - 值越小: 更敏感，检测更多变化点
        /// - 值越大: 更保守，只检测显著变化
        /// 公式: threshold = 3.0 * sensitivityFactor
        /// 注: 基准阈值3.0对应F分布在α=0.05时的临界值</param>
        /// <returns>检测到的变化点列表</returns>
        public static List<ChangePoint> Detect(float[] values, float sensitivityFactor = 2.25f)
        {
            var changePoints = new List<ChangePoint>();

            // 至少需要6个点（前后各3个）才能进行有效的回归分析
            if (values.Length < 6)
                return changePoints;

            // 全数据回归（用于Chow Test）
            float[] fullX = Enumerable.Range(0, values.Length).Select(x => (float)x).ToArray();
            var fullReg = SimpleLinearRegression(fullX, values);
            float rssPooled = fullReg.Residuals.Sum(r => r * r); // 全数据的残差平方和

            // 从第3个点到倒数第3个点遍历，确保两侧都有足够的数据
            for (int i = 3; i < values.Length - 3; i++)
            {
                // 左段：从开始到变化点之前的所有数据
                var leftData = values.Take(i).ToArray();
                float[] leftX = Enumerable.Range(0, leftData.Length).Select(x => (float)x).ToArray();

                // 右段：从变化点到结束的所有数据
                var rightData = values.Skip(i).ToArray();
                float[] rightX = Enumerable.Range(0, rightData.Length).Select(x => (float)x).ToArray();

                // 对左段进行线性回归
                var leftReg = SimpleLinearRegression(leftX, leftData);

                // 对右段进行线性回归
                var rightReg = SimpleLinearRegression(rightX, rightData);

                // Chow Test: 计算F统计量
                // F = [(RSS_pooled - (RSS_left + RSS_right)) / k] / [(RSS_left + RSS_right) / (n - 2k)]
                // 其中 k = 2 (截距 + 斜率), n = 总样本量

                float rssLeft = leftReg.Residuals.Sum(r => r * r);
                float rssRight = rightReg.Residuals.Sum(r => r * r);
                float rssSeparate = rssLeft + rssRight;

                int n = values.Length;
                int k = 2; // 参数数量：截距 + 斜率

                // 避免除零
                if (rssSeparate < MathHelper.MIN_STD_DEV)
                {
                    // 如果分段回归的残差接近0，说明完美拟合，这是异常情况
                    continue;
                }

                // 计算F统计量
                float numerator = (rssPooled - rssSeparate) / k;
                float denominator = rssSeparate / (n - 2 * k);

                if (denominator < MathHelper.MIN_STD_DEV)
                {
                    continue;
                }

                float fStat = numerator / denominator;

                // F统计量必须为正值才有意义（表示分段拟合显著优于全数据拟合）
                if (fStat <= 0)
                    continue;

                // 计算斜率差异（用于判断变化类型和幅度）
                float slopeDiff = Math.Abs(leftReg.Slope - rightReg.Slope);

                // 使用F统计量作为score
                // F值越大，结构突变越显著
                // 典型的F临界值：F(2, n-4, α=0.05) ≈ 3.0 - 4.0
                float score = fStat;

                // 应用灵敏度因子进行检测
                // 基准阈值：F = 3.0 (对应 p ≈ 0.05 的临界值)
                float baseThreshold = 3.0f;
                float threshold = baseThreshold * sensitivityFactor;

                if (score > threshold)
                {
                    // 判断趋势变化方向：右段斜率更大表示趋势加速向上
                    TrendDirection changeType = rightReg.Slope > leftReg.Slope
                        ? TrendDirection.Increasing
                        : TrendDirection.Decreasing;

                    // 计算置信度：使用标准统计方法（基于F分布）
                    // 1. 统计显著性：F统计量相对于阈值的比率
                    float statisticalSignificance = score / threshold;

                    // 2. 效应大小：标准化的斜率差异
                    // 使用两段回归的合并标准误差来标准化斜率差异
                    float pooledSE = (float)Math.Sqrt((rssLeft + rssRight) / (n - 4));

                    // Cohen's d for slope difference
                    float effectSize = slopeDiff / (pooledSE + MathHelper.MIN_STD_DEV);

                    // 归一化效应大小：0.2=小, 0.5=中, 0.8=大
                    effectSize = Math.Min(effectSize / 0.8f, 1.0f);

                    // 3. 综合置信度：使用几何平均数
                    // 确保斜率变化既要统计显著，又要有实际意义
                    float rawConfidence = (float)Math.Sqrt(statisticalSignificance * effectSize);

                    // 4. 映射到[0, MAX_CONFIDENCE]区间
                    float confidence = Math.Min(MAX_CONFIDENCE, 1.0f - (float)Math.Exp(-rawConfidence));

                    changePoints.Add(new ChangePoint(
                        index: i,
                        score: score,
                        changeType: changeType,
                        magnitude: slopeDiff,
                        confidence: confidence
                    ));
                }
            }

            return changePoints;
        }

        /// <summary>
        /// 执行简单线性回归分析
        /// </summary>
        private static (float Slope, float Intercept, float[] FittedValues, float[] Residuals)
            SimpleLinearRegression(float[] x, float[] y)
        {
            int n = x.Length;
            float xMean = x.Average();
            float yMean = y.Average();

            float numerator = 0f;
            float denominator = 0f;

            for (int i = 0; i < n; i++)
            {
                numerator += (x[i] - xMean) * (y[i] - yMean);
                denominator += (x[i] - xMean) * (x[i] - xMean);
            }

            float slope = denominator == 0 ? 0 : numerator / denominator;
            float intercept = yMean - slope * xMean;

            float[] fittedValues = new float[n];
            float[] residuals = new float[n];

            for (int i = 0; i < n; i++)
            {
                fittedValues[i] = slope * x[i] + intercept;
                residuals[i] = y[i] - fittedValues[i];
            }

            return (slope, intercept, fittedValues, residuals);
        }
    }
}
